function [ Fgrid ] = Matrix_Bivariate_4pointRefinement_boundary_V2( grid , omega )
%-----------------------------------------------------------------
% Input:
% grid(m,n) - a grid over Z^2
% Output:
% Fgrid(2*m-1,2*n-1)
%
% ASSUMPTION: m,n > 4
%-----------------------------------------------------------------
% ABSTRACT
% Tensor product, includes boundaries. the middle point is insert in the order y->x of
% calculation
%-----------------------------------------------------------------
% NIr Sharon, 29-05-12
%-----------------------------------------------------------------

m = size(grid,1);
n = size(grid,2);
N = size(grid,3);
Z = zeros(N);
I = eye(N);

Fgrid = zeros(2*m-1,2*n-1,N,N);


%-----------------------------------------------------------
% First main loops: For vertices (new and old) on the grid

for x=1:(m)
    for y=1:(n)
        % Interpolation
        Fgrid(2*x-1,2*y-1,:,:) = squeeze(grid(x,y,:,:));
        % Y direction
        if (y~=n)
            if (y==1)
                Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)) , squeeze(grid(x,y+3,:,:)),omega);
            else
                if (y==(n-1))
                    Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y-2,:,:)),omega);
                else
                    % Divide into cases of zero matrix
                    % left side zero
                    if (isequal(Z,squeeze(grid(x,y-1,:,:))))
                        if ((y+3)<=n)
                            Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)) , squeeze(grid(x,y+3,:,:)),omega);
                        else
                            Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)) , I ,omega);
                        end
                        % right side zero
                    elseif (isequal(Z,squeeze(grid(x,y+2,:,:))))
                        if ((y-2)>=1)
                            Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y-2,:,:)),omega);
                        else
                            Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y-1,:,:)) , I ,omega);
                        end
                        % the main possibility
                    else
                        Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage(squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)),omega);
                    end
                end
            end
        end
        % X direction
        if (x~=m)
            if (x==1)
                Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)) , squeeze(grid(x+3,y,:,:)),omega);
            else
                if (x==(m-1))
                    Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x+1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x-1,y,:,:)) , squeeze(grid(x-2,y,:,:)),omega);
                else
                    % Divide into cases of zero matrix
                    % left side zero
                    if (isequal(Z,squeeze(grid(x-1,y,:,:))))
                        if ((x+3)<=m)
                            Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)) , squeeze(grid(x+3,y,:,:)),omega);
                        else
                            Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)) , I, omega);
                        end
                        % right side zero
                    elseif (isequal(Z,squeeze(grid(x+2,y,:,:))))
                        if ((x-2)>=1)
                            Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x+1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x-1,y,:,:)) , squeeze(grid(x-2,y,:,:)),omega);
                        else
                            Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary_V2(squeeze(grid(x+1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x-1,y,:,:)) , I,omega);
                        end
                        % the main possibility
                    else
                        Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)),omega);
                    end
                end
            end
        end
    end
    
    %-----------------------------------------------------------
    % Second main loops : grid face centers
    
    for x=1:(m-1)
        for y=1:(n-1)
            % Y direction
            if (y==1)
                m_y = Matrix_FourPointAvarage_boundary_V2(squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)) , squeeze(Fgrid(2*x,2*y+5,:,:)),omega);
            else
                if (y==(n-1))
                    m_y = Matrix_FourPointAvarage_boundary_V2(squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y-3,:,:)) , squeeze(Fgrid(2*x,2*y-5,:,:)),omega);
                else
                    % Divide into cases of zero matrix
                    % left side zero
                    if (isequal(Z,squeeze(Fgrid(2*x,2*y-3,:,:))))
                        if ((2*y+5)<=(2*n-1))
                            m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)) , squeeze(Fgrid(2*x,2*y+5,:,:)),omega);
                        else
                            m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)) , I,omega);
                        end
                        % right side zero
                    elseif (isequal(Z,squeeze(Fgrid(2*x,2*y+3,:,:))))
                        if ((2*y-5)>=1)
                            m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y-3,:,:)) , squeeze(Fgrid(2*x,2*y-5,:,:)),omega);
                        else
                            m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y-3,:,:)) , I,omega);
                        end
                        % the main possibility
                    else
                        m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y-3,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)),omega);
                    end
                end
            end
            % X direction
            if (x==1)
                m_x = Matrix_FourPointAvarage_boundary_V2(squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)) , squeeze(Fgrid(2*x+5,2*y,:,:)),omega);
            else
                if (x==(m-1))
                    m_x = Matrix_FourPointAvarage_boundary_V2(squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x-3,2*y,:,:)) , squeeze(Fgrid(2*x-5,2*y,:,:)),omega);
                else
                    % Divide into cases of zero matrix
                    % left side zero
                    if (isequal(Z,squeeze(Fgrid(2*x-3,2*y,:,:))))
                        if ((2*x+5)<=(2*m-1))
                            m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)) , squeeze(Fgrid(2*x+5,2*y,:,:)),omega);
                        else
                            m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)) , I ,omega);
                        end
                        % right side zero
                    elseif (isequal(Z,squeeze(Fgrid(2*x+3,2*y,:,:))))
                        if ((2*x-5)>=1)
                            m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x-3,2*y,:,:)) , squeeze(Fgrid(2*x-5,2*y,:,:)),omega);
                        else
                            m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x-3,2*y,:,:)) , I ,omega);
                        end
                        % the main possibility
                    else
                        m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-3,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)),omega);
                    end
                end
            end
                % dif(x,y) = m_x-m_y;
                Fgrid(2*x,2*y,:,:) = GM2(m_x,m_y,.5);
            end
        end
        
    end
    
    
    
